MAREANO - NiN - hiDens >1031.8

Packages

library(dplyr)
library(ggcorrplot)
library(ggpubr)
library(patchwork)
# library("here")
# library("bookdown")
# library("downloadthis")
library(vegan)
library(plyr)
library(e1071)
library(tidyverse)
library(viridis)
library(GGally)
library(ggrepel)
library(ggordiplots)
library(readr)
library(RColorBrewer)
library(oce) 
library(plotly)


source("gg_ordisurf_viridis.R")

Data

env_sort <- read.csv(file.path(dataPath,"inputs/HDensDeep_env_sort_2022-09-23.csv")) %>% as.data.frame
otu_sort <- read.csv(file.path(dataPath,"inputs/HDensDeep_otu_sort_2022-09-23.csv")) %>% as.data.frame


table(env_sort$SampID2 == otu_sort$SampID)

< table of extent 0 >

#must be same length and equal to unique sample number length - the ordered lists are assumed to be directly relatable on a row by row basis

dim(env_sort)

[1] 117 142

dim(otu_sort)

[1] 117 351

remove any variables with not enough coverage

env_sort <- env_sort %>% select(-c("BO22_lightbotmean_bdmean",
                               "BO22_lightbotltmax_bdmean",
                               "BO22_lightbotltmin_bdmean",
                               "BO22_lightbotrange_bdmean"))

Plot map to check location of samples

env_sort_locations<- ggplot(data = env_sort,
                      aes(x = X.y,
                          y = Y)) +
  theme_classic() +
  geom_point(aes(colour = bathy),
             size = 2) +
  scale_colour_gradient2(low = "red",
                         mid = "yellow",
                         high = "green") +
  ggtitle("Location of samples in sorted env file")

env_sort_locations

Data cleaning

## Removing NAs
otuCompl <- otu_sort[complete.cases(env_sort[, -c(1:23)]), ]
envCompl <- env_sort[complete.cases(env_sort[, -c(1:23)]), ]

## Removing observations with less than 4 OTUs
sel <- rowSums(otuCompl[, -c(1:2)]) >= 4
otuSel <- otuCompl[sel, ]
envSel <- envCompl[sel, ]


## Removing phosphate
#envSel <- envSel %>% select(-phosphate_mean.tif)


dim(otuSel); dim(envSel)

[1]  95 351
[1]  95 138

Show what samples are left after complete cases and >4 OTUs filters

env_sel_locations<- ggplot(data = envSel,
                      aes(x = X.y,
                          y = Y)) +
  theme_classic() +
  geom_point(aes(colour = bathy),
             size = 2) +
  scale_colour_gradient2(low = "red",
                         mid = "yellow",
                         high = "green") +
  ggtitle("Location of samples in sorted env file")

env_sel_locations

Show what got removed in complete cases filter

envCompl_ccrem<-env_sort%>%filter(!SampID%in%envCompl$SampID)

envCompl_ccrem_locations<- ggplot(data = envCompl_ccrem,
                      aes(x = X.y,
                          y = Y)) +
  theme_classic() +
  geom_point(aes(colour = bathy),
             size = 2) +
  scale_colour_gradient2(low = "red",
                         mid = "yellow",
                         high = "green") +
  ggtitle("Location of removed complete case samples resulting in envSel file")

envCompl_ccrem_locations

Show what got removed in <4 OTUs filter

inv.sel <- rowSums(otuCompl[, -c(1:2)]) < 4
env.invSel <- envCompl[inv.sel, ]

env.invSel_locations<- ggplot(data = env.invSel,
                      aes(x = X.y,
                          y = Y)) +
  theme_classic() +
  geom_point(aes(colour = bathy),
             size = 2) +
  scale_colour_gradient2(low = "red",
                         mid = "yellow",
                         high = "green") +
  ggtitle("Location of removed samples with <4 OTUs resulting in envSel file")

env.invSel_locations

[Optional] Data thinning

# otu_red <- otu1[1:500, ]
# otu <- otu_red
# 
# env_red <- env[1:500, ]
# env <- env_red

Whole dataset for first run

otu<-otuSel[,-c(1)]
env<-envSel[, -c(1,3:23)]

table(is.na(otu))

FALSE 
33250 

table(is.na(env))

FALSE 
11020 

#str(otuSel)
#str(envSel)

Splitting in subsets

# ## Class 1
# otu1 <- subset(otuSel, envSel$SplitRev == 1)
# env1 <- envSel %>% filter(SplitRev == 1)
# 
# ## Class 2
# otu2 <- subset(otuSel, envSel$SplitRev == 2)
# env2 <- envSel %>% filter(SplitRev == 2)
# 
# ## Class 4
# otu3 <- subset(otuSel, envSel$SplitRev == 3)
# env3 <- envSel %>% filter(SplitRev == 3)
# 
# ## Class 4
# otu4 <- subset(otuSel, envSel$SplitRev == 4)
# env4 <- envSel %>% filter(SplitRev == 4)
# 
# ## Class 6
# otu6 <- subset(otuSel, envSel$SplitRev == 6)
# env6 <- envSel %>% filter(SplitRev == 6)
# 
# ## Class 7
# otu7 <- subset(otuSel, envSel$SplitRev == 7)
# env7 <- envSel %>% filter(SplitRev == 7)
# 
# ## Class 8
# otu8 <- subset(otuSel, envSel$SplitRev == 8)
# env8 <- envSel %>% filter(SplitRev == 8)

Selecting subset

# ## Selecting 
# otu <- otu2
# env <- env2

format data correctly

#env$coords.x1<-as.numeric(env$coords.x1)
#env$coords.x2<-as.numeric(env$coords.x2)

Abundance weighting

Make function You might have to drop variables that have been imported as character

otu_pa <- decostand(x = otu[, -c(1)],
                    method = "pa")

# y = ax^w         # power transformation formula
dt <- otu[, -c(1:2)]          # species data to transform
x_mn <- min(dt[dt > 0])
x_mx <- max(dt)
rng <- 6           # abundance range
w <- log(rng) / (log(x_mx) - log(x_mn))
a <- x_mn^(-w)
# otu_6 <- a * dt[, -c(1:3)]^w
otu_6 <- a * dt^w
range(otu_6)

[1] 0 6

Ordination methods

DCA PA

Sys.time()

[1] "2022-09-26 11:11:18 CEST"

dca_pa <- decorana(veg = otu_pa)

some species were removed because they were missing in the data

print(dca_pa, head=T)

Call:
decorana(veg = otu_pa) 

Detrended correspondence analysis with 26 segments.
Rescaling of axes with 4 iterations.

                  DCA1   DCA2   DCA3   DCA4
Eigenvalues     0.4911 0.3971 0.3139 0.2894
Decorana values 0.5106 0.4359 0.3537 0.3344
Axis lengths    4.9733 5.3038 3.6768 4.4672

GNMDS PA

Distances - don’t run if loading saved object

## Bray-Curtis
dist_pa <- vegdist(x = otu_pa, method = "bray")

## Geodist
ep <- 0.8     # epsilon
geodist_pa <- isomapdist(dist = dist_pa, epsilon = ep)

Save the result - don’t run if loading saved object

saveRDS(geodist_pa,
        file = (file.path(dataPath,"inputs/HiDens_geodist_pa.rds")))

Ordination

monoMDS - don’t run if loading saved object

Sys.time()

[1] "2022-09-26 11:11:28 CEST"

d <- 2
mds_pa <- list()

for (i in 1:100) {
  mds_pa[[i]]<-monoMDS(geodist_pa,
                    matrix(c(runif(dim(otu_pa)[1]*d)),
                           nrow = dim(otu_pa)[1]),
                    k = d,
                    model = "global",
                    maxit = 2000,
                    smin = 1e-7,
                    sfgrmin = 1e-7)
}

Save the result - don’t run if loading saved object

saveRDS(mds_pa,
        file = file.path(dataPath,"inputs/HiDens_mds_pa.rds"))

Best nmds solution - PA

Make function?

# Loading geodist object
# geodist_nfi <- readRDS(file = "../SplitRev2_geodist_pa_full.rds")

# Loading mds results
# mds <- readRDS(file = "../SplitRev2_mds_pa_full.rds")

## Extracting the stress of each nmds iteration
mds_stress_pa<-unlist(lapply(mds_pa, function(v){v[[22]]})) 

ordered_pa <-order(mds_stress_pa)

## Best, second best, and worst solution
mds_stress_pa[ordered_pa[1]]

[1] 0.2475837

mds_stress_pa[ordered_pa[2]]

[1] 0.2475837

mds_stress_pa[ordered_pa[10]]

[1] 0.2475837

## Scaling of axes to half change units and varimax rotation by postMDS
mds_best_pa<-postMDS(mds_pa[[ordered_pa[1]]],
                  geodist_pa, 
                  pc = TRUE, 
                  halfchange = TRUE, 
                  threshold = ep)     # Is this threshold related to the epsilon above?
mds_best_pa

Call:
monoMDS(dist = geodist_pa, y = matrix(c(runif(dim(otu_pa)[1] *      d)), nrow = dim(otu_pa)[1]), k = d, model = "global", maxit = 2000,      smin = 1e-07, sfgrmin = 1e-07) 

Non-metric Multidimensional Scaling

95 points, dissimilarity ‘bray shortest isomap’, call ‘isomapdist(dist = dist_pa, epsilon = ep)’

Dimensions: 2 
Stress:     0.2475837 
Stress type 1, weak ties
Scores scaled to unit root mean square, rotated to principal components
Stopped after 119 iterations: Stress nearly unchanged (ratio > sratmax)

mds_secbest_pa <- postMDS(mds_pa[[ordered_pa[2]]],
                          geodist_pa, 
                          pc = TRUE, 
                          halfchange = TRUE, 
                          threshold = ep)
mds_secbest_pa

Call:
monoMDS(dist = geodist_pa, y = matrix(c(runif(dim(otu_pa)[1] *      d)), nrow = dim(otu_pa)[1]), k = d, model = "global", maxit = 2000,      smin = 1e-07, sfgrmin = 1e-07) 

Non-metric Multidimensional Scaling

95 points, dissimilarity ‘bray shortest isomap’, call ‘isomapdist(dist = dist_pa, epsilon = ep)’

Dimensions: 2 
Stress:     0.2475837 
Stress type 1, weak ties
Scores scaled to unit root mean square, rotated to principal components
Stopped after 90 iterations: Scale factor of gradient nearly zero (< sfgrmin)

## Procrustes comparisons
procr_pa <- procrustes(mds_best_pa,
                    mds_secbest_pa,
                    permutations=999)
protest(mds_best_pa,
        mds_secbest_pa,
        permutations=999)

Call:
protest(X = mds_best_pa, Y = mds_secbest_pa, permutations = 999) 

Procrustes Sum of Squares (m12 squared):        8.014e-09 
Correlation in a symmetric Procrustes rotation:     1 
Significance:  0.001 

Permutation: free
Number of permutations: 999

plot(procr_pa)

png(file=file.path(dataPath,"outputs/HiDens_procrustes_pa.png"), width=1000, height=700)

plot(procr_pa)
dev.off()

png 
  2 

Correlation of axis: DCA vs NMDS - PA

# Extracting ordination axis
ax <- 2
axis_pa <- cbind(mds_best_pa$points,
                 scores(dca_pa,
                        display = "sites",
                        origin = TRUE)[, 1:ax])

ggcorr(axis_pa, 
       method=c("everything","kendall"), 
       label = TRUE,
       label_size = 3, 
       label_color = "black",  
       nbreaks = 8,
       label_round = 3,
       low = "red",
       mid = "white",
       high = "green")

NA
NA

Save the figure

ggsave(filename = file.path(dataPath,"outputs/HiDens_correlationPCAvsNMDS_PA.png"),
       device = "png",
       dpi=300 )

Saving 7 x 7 in image

# Switching direction of NMDS1
# mds_best$points[, 1] <- -mds_best$points[, 1]

DCA R6

dca_r6 <- decorana(veg = otu_6)

some species were removed because they were missing in the data

print(dca_r6, head=T)

Call:
decorana(veg = otu_6) 

Detrended correspondence analysis with 26 segments.
Rescaling of axes with 4 iterations.

                  DCA1   DCA2   DCA3   DCA4
Eigenvalues     0.5090 0.4363 0.2505 0.2718
Decorana values 0.5278 0.4302 0.3405 0.2632
Axis lengths    4.3041 5.0978 3.5602 2.8604

GNMDS R6

Distances - don’t run if loading saved object

## Bray-Curtis
dist_r6 <- vegdist(x = otu_6, method = "bray")

## Geodist
ep <- 0.80     # epsilon
geodist_r6 <- isomapdist(dist = dist_r6, epsilon = ep)

Save the result - don’t run if loading saved object

saveRDS(geodist_r6,
        file = (file.path(dataPath,"inputs/HiDens_geodist_r6.rds")))

Ordination

200 reps - don’t run if loading saved object

took 20mins

# monoMDS
d <- 2
mds_r6 <- list()

Sys.time()

[1] "2022-09-26 11:12:45 CEST"

for (i in 1:200) {
  mds_r6[[i]]<-monoMDS(geodist_r6,
                    matrix(c(runif(dim(otu_6)[1]*d)),
                           nrow = dim(otu_6)[1]),
                    k = d,
                    model = "global",
                    maxit = 2000,
                    smin = 1e-7,
                    sfgrmin = 1e-7)
}
Sys.time()

[1] "2022-09-26 11:12:48 CEST"

Save the result - don’t run if loading saved object

saveRDS(mds_r6,
        file = file.path(dataPath,"inputs/HiDens_mds_r6.rds")) 

Best nmds solution r6 200 rep

# Loading geodist object
# geodist_nfi <- readRDS(file = "../SplitRev2_geodist_nfi.rds")

# Loading mds results
# mds <- readRDS(file = "../SplitRev2_mds.rds")

## Extracting the stress of each nmds iteration
mds_stress_r6<-unlist(lapply(mds_r6, function(v){v[[22]]})) 

ordered_r6 <-order(mds_stress_r6)

## Best, second best, and worst solution
mds_stress_r6[ordered_r6[1]]

[1] 0.2353797

mds_stress_r6[ordered_r6[2]]

[1] 0.2353799

mds_stress_r6[ordered_r6[100]]

[1] 0.237854

## Scaling of axes to half change units and varimax rotation by postMDS
mds_best_r6<-postMDS(mds_r6[[ordered_r6[1]]],
                     geodist_r6, 
                     pc = TRUE, 
                     halfchange = TRUE, 
                     threshold = ep)     # Is this threshold related to the epsilon above?
mds_best_r6

Call:
monoMDS(dist = geodist_r6, y = matrix(c(runif(dim(otu_6)[1] *      d)), nrow = dim(otu_6)[1]), k = d, model = "global", maxit = 2000,      smin = 1e-07, sfgrmin = 1e-07) 

Non-metric Multidimensional Scaling

95 points, dissimilarity ‘bray shortest isomap’, call ‘isomapdist(dist = dist_r6, epsilon = ep)’

Dimensions: 2 
Stress:     0.2353797 
Stress type 1, weak ties
Scores scaled to unit root mean square, rotated to principal components
Stopped after 103 iterations: Stress nearly unchanged (ratio > sratmax)

mds_secbest_r6<-postMDS(mds_r6[[ordered_r6[2]]],
                        geodist_r6, 
                        pc = TRUE, 
                        halfchange = TRUE, 
                        threshold = ep)
mds_secbest_r6

Call:
monoMDS(dist = geodist_r6, y = matrix(c(runif(dim(otu_6)[1] *      d)), nrow = dim(otu_6)[1]), k = d, model = "global", maxit = 2000,      smin = 1e-07, sfgrmin = 1e-07) 

Non-metric Multidimensional Scaling

95 points, dissimilarity ‘bray shortest isomap’, call ‘isomapdist(dist = dist_r6, epsilon = ep)’

Dimensions: 2 
Stress:     0.2353799 
Stress type 1, weak ties
Scores scaled to unit root mean square, rotated to principal components
Stopped after 69 iterations: Stress nearly unchanged (ratio > sratmax)

## Procrustes comparisons
procr_r6 <- procrustes(mds_best_r6,
                       mds_secbest_r6,
                       permutations=999)
protest(mds_best_r6,
        mds_secbest_r6,
        permutations=999)

Call:
protest(X = mds_best_r6, Y = mds_secbest_r6, permutations = 999) 

Procrustes Sum of Squares (m12 squared):        6.903e-06 
Correlation in a symmetric Procrustes rotation:     1 
Significance:  0.001 

Permutation: free
Number of permutations: 999

plot(procr_r6)

png(file.path(dataPath,"outputs/HiDens_procrustes_r6.png"), width=1000, height=700,) #added 1000

plot(procr_r6)
dev.off()

png 
  2 

#### 1000 reps - don’t run if loading saved object Currently commented out as it gained nothing but took extra time. Can be removed in due course.

Correlation of axis: DCA vs NMDS

retain the 200 rep version

# Extracting ordination axis
ax <- 2
axis_r6 <- cbind(mds_best_r6$points,
                 scores(dca_r6,
                        display = "sites",
                        origin = TRUE)[, 1:ax])

ggcorr(axis_r6, 
       method=c("everything","kendall"), 
       label = TRUE,
       label_size = 3, 
       label_color = "black",  
       nbreaks = 8,
       label_round = 3,
       low = "red",
       mid = "white",
       high = "green")

Save the figure

ggsave(filename = file.path(dataPath,"outputs/HiDens_correlationPCAvsNMDS_r6.png"),
       device = "png",
       dpi=300 )

Saving 7 x 7 in image

# Switching direction of NMDS1
# mds_best$points[, 1] <- -mds_best$points[, 1]

Plotting DCA & GNMDS

## Adding scores to data frame
otu_6$gnmds1 <- mds_best_r6$points[, 1]
otu_6$gnmds2 <- mds_best_r6$points[, 2]
otu_6$dca1 <- scores(dca_r6, display = "sites", origin = TRUE)[, 1]
otu_6$dca2 <- scores(dca_r6, display = "sites", origin = TRUE)[, 2]

p_gnmds_r6 <- ggplot(data = otu_6,
                     aes(x = gnmds1,
                         y = gnmds2)) +
  theme_classic() +
  coord_fixed() +
  ggtitle("GNMDS",
          subtitle = "First run") +
  geom_point(colour = "red") +
  geom_vline(xintercept = 0,
             linetype = 2,
             colour = "lightgray") +
  geom_hline(yintercept = 0,
             linetype = 2,
             colour = "lightgray")

p_dca_r6 <- ggplot(data = otu_6,
                   aes(x = dca1,
                       y = dca2)) +
  theme_classic() +
  coord_fixed() +
  ggtitle("DCA",
          subtitle = "First run") +
  geom_point(colour = "red") +
  geom_vline(xintercept = 0,
             linetype = 2,
             colour = "lightgray") +
  geom_hline(yintercept = 0,
             linetype = 2,
             colour = "lightgray")

p_gnmds_r6 + p_dca_r6

NA
NA

Save the figure

ggsave(file.path(dataPath,"outputs/HiDens_gnmds_dca_r6.png"),
       device = "png", 
       dpi=300)

Saving 7 x 7 in image

Species-environment relationships

Selecting ordination

ord <- mds_best_r6

## Axis scores if selected ord is GNMDS
axis <- ord$points %>% as.data.frame

## Axis scores if selected ord is DCA
# axis <- scores(ord,
#                display = "sites",
#                origin = TRUE)[, 1:ax])

Create additional variables

decided to make MLD-bathy vars

env<-env %>% 
  mutate ("MLDmean_bathy"=MLDmean_Robinson-(bathy*-1),
          "MLDmin_bathy"=MLDmin_Robinson-(bathy*-1),
          "MLDmax_bathy"=MLDmax_Robinson-(bathy*-1))

env$MLDmean_bathy<-cut(env$MLDmean_bathy, 
      breaks=c(-2560, -20,20,130),#checked range of values first (min -2554, max 123)
      labels=c('belowMLD','onPycno','inMixLayer'))
env$MLDmin_bathy<-cut(env$MLDmin_bathy, 
      breaks=c(-2560, -20,20,130),#checked range of values first (min -2554, max 123)
      labels=c('belowMLD','onPycno','inMixLayer'))
env$MLDmax_bathy<-cut(env$MLDmax_bathy, 
      breaks=c(-2560, -20,20,130),#checked range of values first (min -2554, max 123)
      labels=c('belowMLD','onPycno','inMixLayer'))

env$swDensRob_avs<-swRho(salinity=env$Smean_Robinson,
                         temperature=env$Tmean_Robinson,
                         pressure=(env$bathy*-1),
                         eos="unesco")

Correlation ordination axes and environmental variables

Removing non-env vars

env_cont<-env%>% select(-c(landscape,sedclass,gmorph, MLDmean_bathy, MLDmax_bathy, MLDmin_bathy, #categorical
                           optional, #not a var
                           MLDmax_Robinson, MLDmean_Robinson, MLDmin_Robinson,  #replaced by new vars
                           MLDsd_Robinson #not meaningful
                           ))
env_cont<-env_cont%>% mutate_if(is.integer,as.numeric)
env_corr <- env_cont %>% select(-c(SampID))

# env_corr$coords.x1<-as.numeric(env_corr$coords.x1)
# env_corr$coords.x2<-as.numeric(env_corr$coords.x2)

env_corr[(!is.numeric(env_corr)),]

Correlations

# Vector to hold correlations
cor_ax1 <- NULL
cor_ax2 <- NULL
pv_ax1 <- NULL
pv_ax2 <- NULL

# NMDS1
for( i in seq(length(env_corr))) {
  ct.i <- cor.test(axis$MDS1,
                   env_corr[, i],
                   method = "kendall")
  cor_ax1[i] <- ct.i$estimate
  pv_ax1[i] <- ct.i$p.value
}

the standard deviation is zerothe standard deviation is zero

# NMDS2
for( i in seq(length(env_corr))) {
  ct.i <- cor.test(axis$MDS2,
                   env_corr[, i],
                   method = "kendall")
  cor_ax2[i] <- ct.i$estimate
  pv_ax2[i] <- ct.i$p.value
}

the standard deviation is zerothe standard deviation is zero

cor_tab <- data.frame(env = names(env_corr),
                      ord_ax1 = cor_ax1,
                      pval_ax1 = pv_ax1,
                      ord_ax2 = cor_ax2,
                      pval_ax2 = pv_ax2)

cor_tab

write.csv(x = cor_tab,
          file = file.path(dataPath,"inputs/HiDens_cor-table_r6_200rep_MLD-bathy.csv"),
          row.names = FALSE)

Dot chart to check for gaps in correlation

cor_a1_sort<-cor_tab%>%
  mutate(abs_ord_ax1=abs(ord_ax1),
         abs_ord_ax2=abs(ord_ax2)) %>%
  arrange(desc(abs_ord_ax1))

cor_a2_sort<-cor_tab%>%
  mutate(abs_ord_ax1=abs(ord_ax1),
         abs_ord_ax2=abs(ord_ax2)) %>%
  arrange(desc(abs_ord_ax2))

dotchart(cor_a1_sort$abs_ord_ax1, main="Absolute (+/-) correlations between envVars and gnmds axis 1")

cor_cut<-0.42 #decide

cor_sel<-subset(cor_a1_sort,abs_ord_ax1>cor_cut)
cor_sel

as.data.frame(cor_sel$env)

Sel env var (top corr) ax1

env_os <- env[, cor_sel$env]
env_os

str(env_os)

'data.frame':   95 obs. of  15 variables:
 $ cobB                      : int  0 0 0 20 0 0 0 0 0 0 ...
 $ temp_max                  : num  5.97 6.76 6.77 1.04 1.06 ...
 $ salt_min                  : num  34.9 34.9 34.9 34.9 34.9 ...
 $ BO22_phosphateltmin_bdmean: num  0.96 0.963 0.946 0.944 0.928 ...
 $ BO22_nitrateltmin_bdmean  : num  13.8 13.8 13.6 13.6 13.4 ...
 $ BO22_phosphatemean_bdmean : num  0.979 0.982 0.968 0.961 0.947 ...
 $ BO22_phosphaterange_bdmean: num  0.979 0.982 0.968 0.961 0.947 ...
 $ BO22_nitratemean_bdmean   : num  14 14.1 13.9 13.9 13.7 ...
 $ BO22_nitraterange_bdmean  : num  14 14.1 13.9 13.9 13.7 ...
 $ BO22_nitrateltmax_bdmean  : num  14.3 14.3 14.2 14.1 13.9 ...
 $ BO22_phosphateltmax_bdmean: num  0.997 1 0.989 0.976 0.964 ...
 $ BO22_silicateltmin_bdmean : num  9.4 9.51 9.07 9.15 8.86 ...
 $ BO22_silicatemean_bdmean  : num  9.82 9.95 9.54 9.47 9.21 ...
 $ BO22_silicaterange_bdmean : num  9.82 9.95 9.54 9.47 9.21 ...
 $ BO22_silicateltmax_bdmean : num  10.25 10.37 10.01 9.77 9.53 ...

Ordisurfs top corr ax 1

ordsrfs <- list(length = ncol(env_os))

for (i in seq(ncol(env_os))) {
  os.i <- gg_ordisurf(ord = ord,
                      env.var = env_os[, i],
                      pt.size = 1,
                      # binwidth = 0.05,
                      var.label = names(env_os)[i],
                      gen.text.size = 10,
                      title.text.size = 15,
                      leg.text.size = 10)
  
  ordsrfs[[i]] <- os.i$plot
}

ordsrfs_plt <- ggarrange(plotlist = ordsrfs,
                         nrow = 5,
                         ncol = 3)

ordsrfs_plt

Save some outputs

ggexport(ordsrfs_plt,
          filename = file.path(dataPath,"outputs/HiDens_ordisurfs_top_corr.png"),
          width = 1000,
          height = 1500)

Dot chart to check for gaps in correlation - AXIS 2

dotchart(cor_a2_sort$abs_ord_ax2, main="Absolute (+/-) correlations between envVars and gnmds axis 2")

cor_cut_a2<-0.26 #decide

cor_sel_a2<-subset(cor_a2_sort,abs_ord_ax2>cor_cut_a2)
cor_sel_a2

as.data.frame(cor_sel_a2$env)

Sel env var (top corr) ax1

env_os_a2 <- env[, cor_sel_a2$env]
env_os_a2

str(env_os_a2)

'data.frame':   95 obs. of  10 variables:
 $ gravel                      : num  18 0.5 0.5 15 15 0 15 0.5 0 0 ...
 $ Vmin_Robinson               : num  -0.059872 -0.051505 -0.027461 -0.000579 0 ...
 $ mud                         : num  24 69.5 69.5 35 60 96 60 69.5 0 0 ...
 $ BO22_carbonphytomean_bdmean : num  0.022 0.0251 0.0215 0.0215 0.0225 ...
 $ BO22_carbonphytorange_bdmean: num  0.022 0.0251 0.0215 0.0215 0.0225 ...
 $ Vmean_Robinson              : num  -0.033264 -0.027986 -0.011453 -0.000508 0 ...
 $ BO22_dissoxltmin_bdmean     : num  260 261 261 306 305 ...
 $ BO22_carbonphytoltmax_bdmean: num  0.0239 0.0302 0.0235 0.0239 0.0258 ...
 $ BO22_chloltmin_ss           : num  0.1638 0.1609 0.1604 0.0342 0.0343 ...
 $ BO22_carbonphytoltmin_ss    : num  0.5228 0.4926 0.4864 0.0942 0.0946 ...

Ordisurfs top corr ax 1

ordsrfs_a2 <- list(length = ncol(env_os_a2))

for (i in seq(ncol(env_os_a2))) {
  os.i_a2 <- gg_ordisurf(ord = ord,
                      env.var = env_os_a2[, i],
                      pt.size = 1,
                      # binwidth = 0.05,
                      var.label = names(env_os_a2)[i],
                      gen.text.size = 10,
                      title.text.size = 15,
                      leg.text.size = 10)
  
  ordsrfs_a2[[i]] <- os.i_a2$plot
}

ordsrfs_plt_a2 <- ggarrange(plotlist = ordsrfs_a2,
                         nrow = 4,
                         ncol = 3)

stat_contour(): Zero contours were generatedno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infstat_contour(): Zero contours were generatedno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infstat_contour(): Zero contours were generatedno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infstat_contour(): Zero contours were generatedno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Inf

ordsrfs_plt_a2

Save some outputs

ggexport(ordsrfs_plt_a2,
          filename = file.path(dataPath,"outputs/HiDens_ordisurfs_top_corr_a2.png"),
          width = 1000,
          height = 1000)

Sel env var (manual)

env_os_m <- env[,c("Tmean_Robinson", #top corr
                  "salt_max", #top corr
                  "Smax_Robinson", #comparison to top corr
                  "swDensRob_avs", #top corr
                  "BO22_icecoverltmax_ss",#top corr ax2
                  "BO22_icecovermean_ss",#top corr
                  "BO22_dissoxmean_bdmean",#top corr
                  #"BO22_carbonphytoltmin_bdmean",#top corr - no clear gradient in ordisurf
                  "BO22_ppltmin_ss", #top corr
                  "X.y", #comparison to Y
                  "Y", #top corr
                  "spd_std", #top corr ax2 (blended model)
                  "CSpdsd_Robinson", #comparison to top corr ax2 (blended model)
                  "mud", #highest sed var ax1 + corr
                  "gravel",#highest sed var ax1 - corr
                  "BO22_silicateltmax_bdmean", #just under top corr ax1
                  "bathy", #intuitive for comparisons
                  "cobB" #added as important, compare to mud /grav
                )]
env_os_m

str(env_os_m)

'data.frame':   95 obs. of  17 variables:
 $ Tmean_Robinson           : num  0.847 1.102 1.49 1.918 2.415 ...
 $ salt_max                 : num  35.1 35.1 35.2 35 35 ...
 $ Smax_Robinson            : num  35 35 35 35.1 35.1 ...
 $ swDensRob_avs            : num  1033 1032 1032 1033 1032 ...
 $ BO22_icecoverltmax_ss    : num  0 0 0 0.00107 0.00141 ...
 $ BO22_icecovermean_ss     : num  0 0 0 0.000171 0.000207 ...
 $ BO22_dissoxmean_bdmean   : num  285 285 285 321 321 ...
 $ BO22_ppltmin_ss          : num  0.000464 0.000407 0.000387 0 0 ...
 $ X.y                      : num  3061 13861 23661 397461 401461 ...
 $ Y                        : num  7053934 7063334 7075134 8609334 8599334 ...
 $ spd_std                  : num  0.0922 0.1036 0.0855 0.0864 0.0894 ...
 $ CSpdsd_Robinson          : num  0.02721 0.02128 0.01167 0.00343 0.00228 ...
 $ mud                      : num  24 69.5 69.5 35 60 96 60 69.5 0 0 ...
 $ gravel                   : num  18 0.5 0.5 15 15 0 15 0.5 0 0 ...
 $ BO22_silicateltmax_bdmean: num  10.25 10.37 10.01 9.77 9.53 ...
 $ bathy                    : num  -953 -869 -835 -963 -943 ...
 $ cobB                     : int  0 0 0 20 0 0 0 0 0 0 ...

Ordisurfs top corr

ordsrfs_m <- list(length = ncol(env_os_m))

for (i in seq(ncol(env_os_m))) {
  os.i_m <- gg_ordisurf(ord = ord,
                      env.var = env_os_m[, i],
                      pt.size = 1,
                      # binwidth = 0.05,
                      var.label = names(env_os_m)[i],
                      gen.text.size = 10,
                      title.text.size = 15,
                      leg.text.size = 10)
  
  ordsrfs_m[[i]] <- os.i_m$plot
}

ordsrfs_plt_m <- ggarrange(plotlist = ordsrfs_m,
                         nrow = 6,
                         ncol = 3)

stat_contour(): Zero contours were generatedno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infstat_contour(): Zero contours were generatedno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infstat_contour(): Zero contours were generatedno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Inf

ordsrfs_plt_m

Save some outputs

ggexport(ordsrfs_plt_m,
          filename = file.path(dataPath,"outputs/HiDens_ordisurfs_man_sel.png"),
          width = 1000,
          height = 1500)

Envfit

## Select if any var should be excluded from envfit (makes less busy to read)
env_os_m_envfit<-env_os_m [,c("Tmean_Robinson", #top corr
                  "salt_max", #top corr
                  "Smax_Robinson", #comparison to top corr
                  "swDensRob_avs", #top corr
                 # "BO22_icecoverltmax_ss",#top corr ax2
                  #"BO22_icecovermean_ss",#top corr
                  "BO22_dissoxmean_bdmean",#top corr
                  #"BO22_carbonphytoltmin_bdmean",#top corr - no clear gradient in ordisurf
                  "BO22_ppltmin_ss", #top corr
                 # "X.y", #comparison to Y
                  #"Y", #top corr
                 # "spd_std", #top corr ax2 (blended model)
                  "CSpdsd_Robinson", #comparison to top corr ax2 (blended model)
                  "mud", #highest sed var ax1 + corr
                  "gravel",#highest sed var ax1 - corr
                  "BO22_silicateltmax_bdmean", #just under top corr ax1
                  "bathy", #intuitive for comparisons
                  "cobB" #added as important, compare to mud /grav
                
  
)]

colnames(env_os_m_envfit)<-c("T", #top corr
                  "sMx____", #top corr
                  "SmaxR", #comparison to top corr
                  "__swDensR", #top corr
                 # "icecovmax",#top corr ax2
                  #"icecovav",#top corr
                  "O2",#top corr
                  #"BO22_carbonphytoltmin_bdmean",#top corr - no clear gradient in ordisurf
                  "ppltmin", #top corr
                 # "X", #comparison to Y
                  #"Y", #top corr
                 # "spd_std", #top corr ax2 (blended model)
                  "CSsd", #comparison to top corr ax2 (blended model)
                  "mud", #highest sed var ax1 + corr
                  "gravel",#highest sed var ax1 - corr
                  "SiMx", #just under top corr ax1
                  "bathy", #intuitive for comparisons
                  "cobB" #added as important, compare to mud /grav
                 )

## Envfot plot
gg_envfit(ord = ord,
          env = env_os_m_envfit,
          pt.size = 1)

## Envfit analysis

ef <- envfit(ord = ord,
             env = env_os_m_envfit,
             # na.rm = TRUE
             )
efDF <- as.data.frame(scores(ef,
                             display = "vectors"))

Save the plot

ggsave(filename = file.path(dataPath,"outputs/HiDens_EnvFit_man_sel_cln.png"),
       device = "png",
       dpi=300 )

Saving 7 x 7 in image

Categorical envVar visualised on the mdsplots

use dataset with categorical var included

env_vis<-env
env_vis$gnmds1 <- otu_6$gnmds1
env_vis$gnmds2 <- otu_6$gnmds2

gnmds w mld mean - bathy

p_mld <- ggplot(data = env_vis,
                     aes(x = gnmds1,
                         y = gnmds2)) +
  theme_classic() +
  coord_fixed() +
  ggtitle("GNMDS coloured by proximity to mixed layer depth",
          subtitle = "First run") +
  geom_point(aes(colour = MLDmean_bathy)) +
  geom_vline(xintercept = 0,
             linetype = 2,
             colour = "lightgray") +
  geom_hline(yintercept = 0,
             linetype = 2,
             colour = "lightgray")

p_mld

Save the plot

ggsave(filename = file.path(dataPath,"outputs/HiDens_gnmds_MLDmeanBathy.png"),
       device = "png",
       dpi=300 )

gnmds w sedclass

deal with sedclass codes

env_vis<- env_vis %>% 
  mutate(
    sedclassName = case_when(
            sedclass == "1" ~ "SedCoverR",
            sedclass == "5" ~ "Rock",
            sedclass == "20" ~ "Mud",
            sedclass == "21" ~ "MwBlock",
            sedclass == "40" ~ "sMud",
            sedclass == "80" ~ "mSand",
            sedclass == "100" ~ "Sand",
            sedclass == "110" ~ "gMud",
            sedclass == "115" ~ "gsMud",
            sedclass == "120" ~ "gmSand",
            sedclass == "130" ~ "gSand",
            sedclass == "150" ~ "MSG",
            sedclass == "160" ~ "sGravel",
            sedclass == "170" ~ "Gravel",
            sedclass == "175" ~ "GravBlock",
            sedclass == "185" ~ "SGBmix",
            sedclass == "205" ~ "S/MwB",
            sedclass == "206" ~ "S/MwG/B",
            sedclass == "215" ~ "SGBalt",
            sedclass == "300" ~ "HardSed",
            sedclass == "500" ~ "Biogenic"
                        
    )
  )

colour palette to cope with up to 25 categorical colours

c25 <- c(
  "dodgerblue2", "#E31A1C", # red
  "green4",
  "#6A3D9A", # purple
  "#FF7F00", # orange
  "black", "gold1",
  "skyblue2", "#FB9A99", # lt pink
  "palegreen2",
  "#CAB2D6", # lt purple
  "#FDBF6F", # lt orange
  "gray70", "khaki2",
  "maroon", "orchid1", "deeppink1", "blue1", "steelblue4",
  "darkturquoise", "green1", "yellow4", "yellow3",
  "darkorange4", "brown"
)

p_sed <- ggplot(data = env_vis,
                     aes(x = gnmds1,
                         y = gnmds2)) +
  theme_classic() +
  coord_fixed() +
  ggtitle("GNMDS coloured by sediment class",
          subtitle = "First run") +
  geom_point(aes(colour = factor(sedclassName))) +
   scale_colour_manual(values=c25)+
  geom_vline(xintercept = 0,
             linetype = 2,
             colour = "lightgray") +
  geom_hline(yintercept = 0,
             linetype = 2,
             colour = "lightgray")+
  guides(colour=guide_legend(ncol=2))

p_sed

Save the plot

ggsave(filename = file.path(dataPath,"outputs/HiDens_gnmds_sedclass.png"),
       device = "png",
       dpi=300 )

gnmds w lanscape

deal with sedclass codes

env_vis<- env_vis %>% 
  mutate(
    landscapeName = case_when(
            landscape == "1" ~ "Strandflat",
            landscape == "21" ~ "ContSlope",
            landscape == "22" ~ "Canyon",
            landscape == "31" ~ "Valley",
            landscape == "32" ~ "Fjord",
            landscape == "41" ~ "DeepPlain",
            landscape == "42" ~ "SlopePlain",
            landscape == "43" ~ "ShelfPlain",
            landscape == "431" ~ "shallowValley"
    )
  )

p_land <- ggplot(data = env_vis,
                     aes(x = gnmds1,
                         y = gnmds2)) +
  theme_classic() +
  coord_fixed() +
  ggtitle("GNMDS coloured by landscape class",
          subtitle = "First run") +
  geom_point(aes(colour = factor(landscapeName))) +
  geom_vline(xintercept = 0,
             linetype = 2,
             colour = "lightgray") +
  geom_hline(yintercept = 0,
             linetype = 2,
             colour = "lightgray")+
  guides(colour=guide_legend(ncol=2))

p_land

Save the plot

ggsave(filename = file.path(dataPath,"outputs/HiDens_gnmds_landscape.png"),
       device = "png",
       dpi=300 )

Gnmds w gmorph

p_gmo <- ggplot(data = env_vis,
                     aes(x = gnmds1,
                         y = gnmds2)) +
  theme_classic() +
  coord_fixed() +
  ggtitle("GNMDS coloured by landscape class",
          subtitle = "First run") +
  geom_point(aes(colour = factor(gmorph))) +
  geom_vline(xintercept = 0,
             linetype = 2,
             colour = "lightgray") +
  geom_hline(yintercept = 0,
             linetype = 2,
             colour = "lightgray")+
  guides(colour=guide_legend(ncol=2))

p_gmo

Save the plot

ggsave(filename = file.path(dataPath,"outputs/HiDens_gnmds_gmorph.png"),
       device = "png",
       dpi=300 )

cat_var_plots<-p_mld+p_sed+p_land+p_gmo

Save the plot

ggexport(cat_var_plots,
          filename = file.path(dataPath,"outputs/HiDens_gnmds_catvar.png"),
          width = 1000,
          height = 800)

Sample identification in the mdsplot

p_gmo <- ggplot(data = env_vis,
                     aes(x = gnmds1,
                         y = gnmds2)) +
  theme_classic() +
  coord_fixed() +
  ggtitle("GNMDS coloured by landscape class",
          subtitle = "First run") +
  geom_point(aes(colour = factor(SampID))) +
  geom_vline(xintercept = 0,
             linetype = 2,
             colour = "lightgray") +
  geom_hline(yintercept = 0,
             linetype = 2,
             colour = "lightgray")+
  guides(colour="none", size="none")

ggplotly(p_gmo)